---
title: "R22.RobustSample"
---

# Load Packages
```{r}
library(tidyverse)
library(estimatr)
```
# Load Data
```{r}
load("Data/data_final.Rdata")
load("Data/gnr_panel.Rdata")
load("Data/ideo_list.Rdata")
```

# Set Covariates
```{r}
covs_list <- c("tx_Hi", "as.factor(DistN)", "ln_PopDens_70", "Cities_10km", "ln_Confs_All", "PrRecH_Adm3", "Isolados1950", "PrWheatAr_1972", "LandIneqRat", "TerrainSlopeIndex", "MedAltitude", "Communist_1975", "MinDistPCP_Map", "log_gnr_cumul_23_Z")
```

```{r}
elections <- c("AC75", "AR76", "AR79", "AR80", "AR83", "AR85", "AR87")
years <- data.frame(election = c("AC75", "AR76", "AR79", "AR80", "AR83", "AR85", "AR87"), 
           year = as.factor(c(1975, 1976, 1979, 1980, 1983, 1985, 1987)))
```

# Drop Parishes, based on Cumulative GNR Events up to Election (App. Figures 8bi)
```{r}
gnr_inclusion_list <- list()
gnr_inclusion_list[[1]] <- data_final %>% mutate(included = 1)
gnr_inclusion_list[[2]] <- data_final %>% filter(is.na(drop_by) | drop_by %in% c("AR79", "AR80", "AR83", "AR85", "AR87")) %>% mutate(included = 2) # AR76 vs AC75
gnr_inclusion_list[[3]] <- data_final %>% filter(is.na(drop_by) | drop_by %in% c("AR80", "AR83", "AR85", "AR87")) %>% mutate(included = 3) # AR79 vs AC75
gnr_inclusion_list[[4]] <- data_final %>% filter(is.na(drop_by) | drop_by %in% c("AR83", "AR85", "AR87")) %>% mutate(included = 4) # AR80 vs AC75
gnr_inclusion_list[[5]] <- data_final %>% filter(is.na(drop_by) | drop_by %in% c("AR85", "AR87")) %>% mutate(included = 5) # AR83 vs AC75
gnr_inclusion_list[[6]] <- data_final %>% filter(is.na(drop_by) | drop_by %in% c("AR87")) %>% mutate(included = 6) # AR85 vs AC75
gnr_inclusion_list[[7]] <- data_final %>% filter(is.na(drop_by)) %>% mutate(included = 7) # AR87 vs AC75
```

# Plot Parish Loss Over Time
```{r}
gnr_dropby_coverage <- list()
gnr_dropby_coverage[[1]] <- data_final %>% select(DICOFRE_Z, exp_cat2) %>% mutate(included = 1)
gnr_dropby_coverage[[2]] <- data_final %>% select(DICOFRE_Z, exp_cat2) %>% filter(DICOFRE_Z %in% c(gnr_inclusion_list[[2]]$DICOFRE_Z)) %>% mutate(included = 2)
gnr_dropby_coverage[[3]] <- data_final %>% select(DICOFRE_Z, exp_cat2) %>% filter(DICOFRE_Z %in% c(gnr_inclusion_list[[3]]$DICOFRE_Z)) %>% mutate(included = 3)
gnr_dropby_coverage[[4]] <- data_final %>% select(DICOFRE_Z, exp_cat2) %>% filter(DICOFRE_Z %in% c(gnr_inclusion_list[[4]]$DICOFRE_Z)) %>% mutate(included = 4)
gnr_dropby_coverage[[5]] <- data_final %>% select(DICOFRE_Z, exp_cat2) %>% filter(DICOFRE_Z %in% c(gnr_inclusion_list[[5]]$DICOFRE_Z)) %>% mutate(included = 5)
gnr_dropby_coverage[[6]] <- data_final %>% select(DICOFRE_Z, exp_cat2) %>% filter(DICOFRE_Z %in% c(gnr_inclusion_list[[6]]$DICOFRE_Z)) %>% mutate(included = 6)
gnr_dropby_coverage[[7]] <- data_final %>% select(DICOFRE_Z, exp_cat2) %>% filter(DICOFRE_Z %in% c(gnr_inclusion_list[[7]]$DICOFRE_Z)) %>% mutate(included = 7)

do.call(rbind, gnr_dropby_coverage) %>% 
  left_join(years %>% mutate(included = seq_along(election))) %>% 
  group_by(exp_cat2, election) %>% count %>% 
  ggplot(aes(x=election, y=n, col = as.factor(exp_cat2), group = as.factor(exp_cat2), shape = as.factor(exp_cat2) )) + geom_point() + geom_line() + theme_light() +
  xlab("Election") + ylab("Number of Remaining Parishes") + 
  scale_color_manual(name = "Land Reform Intensity", labels = c("No LR", "LR < 50%, or adjacent", "LR >= 50%"),
                     values = c("#252525", "#636363", "#969696"))+
  scale_shape_manual(name = "Land Reform Intensity", labels = c("No LR", "LR < 50%, or adjacent", "LR >= 50%"),
                       values = c(16, 17, 15)) 

ggsave(plot = last_plot(), 
       "Figures/Afg8b1.1_RobustGNRDrop_Trend.pdf")
```

```{r}

baseline_list <- list()  
for (i in 2:7){ # for elections between 1976 and 1987

aggregated_1 <- ideo_list[[1]] %>% select(-election)   # Get vote shares for baseline year (1975)
colnames(aggregated_1)[3] <- paste(colnames(aggregated_1)[3], "1", sep = "_")
 
aggregated_2 <- ideo_list[[i]] %>% select(-election) # Get vote shares for year of interest
colnames(aggregated_2)[3] <- paste(colnames(aggregated_2)[3], "2", sep = "_")

# First Difference
aggregated_fd <- aggregated_2 %>% inner_join(aggregated_1) %>% rowwise() %>% # Subtract 1975 vote shares from vote shares of year of interest
  mutate(FD = votes_2 - votes_1) %>%
  select(DICOFRE_Z, ideology2, FD) %>%
  pivot_wider(c(DICOFRE_Z), names_from = "ideology2", values_from = "FD")

data_sample_g <- data_final %>% left_join(aggregated_fd) %>% # Join vote share first difference and cumulative GNR events up to election year to main data set
  left_join(gnr_panel %>% filter(elec_number == i)) %>%  
  filter(DICOFRE_Z %in% c(gnr_inclusion_list[[i]]$DICOFRE_Z)) # Drop parishes based on median cumulative GNR visit by election year

output <- list() # Subset for (3) groups (land reform intensity)
for (g in 0:2){ 
 
output[[g+1]] <- data.frame(election = elections[i], group = g,  model = colnames(aggregated_fd)[-c(1)], tx_Hi = NA, se = NA, n = NA) %>% filter(model != "NA")

      # Run models for each ideology category (4)
      for (r in 1:nrow(output[[g+1]])){ # r<-1
      m.1 <- lm_robust(as.formula(paste(output[[g+1]]$model[r], paste(c(covs_list), collapse= " + "), sep = " ~ ")),
                data = data_sample_g %>% filter(exp_cat2 == g), weight = weights, cluster = subclass)
      
           output[[g+1]][r, "tx_Hi"] <- m.1$coefficients["tx_Hi"]
           output[[g+1]][r, "se"] <- m.1$std.error["tx_Hi"]
           output[[g+1]][r, "n"] <- length(m.1$fitted.values)
      }

}

output_list <- do.call(rbind, output)
baseline_list[[i]] <- output_list

}

# Plot Results

baseline_list_all <- do.call(rbind, baseline_list) %>% 
  mutate(model = ifelse(model == "Center", "Center-Left", model), 
         model = ifelse(model == "Left", "Far Left", model))

baseline_list_all <- baseline_list_all %>% left_join(years)
baseline_list_all$model <- factor(baseline_list_all$model, level = c("Communist", "Far Left", "Center-Left", "Right"))

ggplot(baseline_list_all, aes(x = year, col = as.factor(group))) +
  geom_point(aes(y = tx_Hi, shape = as.factor(group)), position = position_dodge(width = 1)) +
  geom_linerange(aes(ymin=tx_Hi - (qt(0.975, n)*se), ymax=tx_Hi + (qt(0.975, n)*se)), position = position_dodge(width = 1), size = 0.5) +
  geom_linerange(aes(ymin=tx_Hi - (qt(0.95, n)*se), ymax=tx_Hi + (qt(0.95, n)*se)), position = position_dodge(width = 1), size = 0.75) +
  geom_hline(aes(yintercept = 0), col = "red") +
  facet_wrap(~model, scales = "free_y") +
  xlab("Election Year") + ylab("Treatment Effect") +
  scale_color_manual(name = "Land Reform Intensity", labels = c("No LR", "LR < 50%, or adjacent", "LR >= 50%"),
                     values = c("#252525", "#636363", "#969696"))+
  scale_shape_manual(name = "Land Reform Intensity", labels = c("No LR", "LR < 50%, or adjacent", "LR >= 50%"),
                       values = c(16, 17, 15)) +
  theme_light() +
  theme(axis.text.x = element_text(angle = 90))

  ggsave(plot = last_plot(),
             "Figures/Afg8b1.2_RobustGNRDrop.pdf",
             width = 9, height = 5)

```

# Drop Parishes that had a Majority of Expripriated Land Returned (App. Figures 8bii)
```{r}
  
baseline_list <- list()  
for (i in 2:7){ # for elections between 1976 and 1987

aggregated_1 <- ideo_list[[1]] %>% select(-election)   # Get vote shares for baseline year (1975)
colnames(aggregated_1)[3] <- paste(colnames(aggregated_1)[3], "1", sep = "_")
 
aggregated_2 <- ideo_list[[i]] %>% select(-election)  # Get vote shares for year of interest
colnames(aggregated_2)[3] <- paste(colnames(aggregated_2)[3], "2", sep = "_")

# First Difference
aggregated_fd <- aggregated_2 %>% inner_join(aggregated_1) %>% rowwise() %>% # Subtract 1975 vote shares from vote shares of year of interest
  mutate(FD = votes_2 - votes_1) %>%
  select(DICOFRE_Z, ideology2, FD) %>%
  pivot_wider(c(DICOFRE_Z), names_from = "ideology2", values_from = "FD")

data_sample_g <- data_final %>% left_join(aggregated_fd) %>% # Join vote share first difference and cumulative GNR events up to election year to main data set
  left_join(gnr_panel %>% filter(elec_number == i)) %>%  
  filter(is.na(pr_haret_exp_coop_85) | pr_haret_exp_coop_85 < 100) # Drop parishes that had no land returned/expropriated or all land returned

output <- list() # Subset for (3) groups (land reform intensity)
for (g in 0:2){ 
 
output[[g+1]] <- data.frame(election = elections[i], group = g,  model = colnames(aggregated_fd)[-c(1)], tx_Hi = NA, se = NA, n = NA) %>% filter(model != "NA")

      # Run models for each ideology category (4)
      for (r in 1:nrow(output[[g+1]])){ 
      m.1 <- lm_robust(as.formula(paste(output[[g+1]]$model[r], paste(c(covs_list), collapse= " + "), sep = " ~ ")),
                data = data_sample_g %>% filter(exp_cat2 == g), weight = weights, cluster = subclass)
      
           output[[g+1]][r, "tx_Hi"] <- m.1$coefficients["tx_Hi"]
           output[[g+1]][r, "se"] <- m.1$std.error["tx_Hi"]
           output[[g+1]][r, "n"] <- length(m.1$fitted.values)
      }

}

output_list <- do.call(rbind, output)
baseline_list[[i]] <- output_list

}

# Plot Results

baseline_list_all <- do.call(rbind, baseline_list) %>% 
  mutate(model = ifelse(model == "Center", "Center-Left", model), 
         model = ifelse(model == "Left", "Far Left", model))

baseline_list_all <- baseline_list_all %>% left_join(years)
baseline_list_all$model <- factor(baseline_list_all$model, level = c("Communist", "Far Left", "Center-Left", "Right"))

ggplot(baseline_list_all, aes(x = year, col = as.factor(group))) +
  geom_point(aes(y = tx_Hi, shape = as.factor(group)), position = position_dodge(width = 1)) +
  geom_linerange(aes(ymin=tx_Hi - (qt(0.975, n)*se), ymax=tx_Hi + (qt(0.975, n)*se)), position = position_dodge(width = 1), size = 0.5) +
  geom_linerange(aes(ymin=tx_Hi - (qt(0.95, n)*se), ymax=tx_Hi + (qt(0.95, n)*se)), position = position_dodge(width = 1), size = 0.75) +
  geom_hline(aes(yintercept = 0), col = "red") +
  facet_wrap(~model, scales = "free_y") +
  xlab("Election Year") + ylab("Treatment Effect") +
  scale_color_manual(name = "Land Reform Intensity", labels = c("No LR", "LR < 50%, or adjacent", "LR >= 50%"),
                     values = c("#252525", "#636363", "#969696"))+
  scale_shape_manual(name = "Land Reform Intensity", labels = c("No LR", "LR < 50%, or adjacent", "LR >= 50%"),
                       values = c(16, 17, 15)) +
  theme_light() +
  theme(axis.text.x = element_text(angle = 90))

  ggsave(plot = last_plot(),
             "Figures/Afg8b2_RobustAllRetDrop.pdf",
             width = 9, height = 5)

```

